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Abstract. We study the dynamical correlation functions and heat conduction for the simplest model of 
quasi one-dimensional (Id) dielectric crystal i.e. a chain of classical particles coupled by quadratic and 
cubic intersite potential. Even in the weakly nonlinear regime, numerical simulation on long enough chains 
reveal sizeable deviations from the perturbative results in the form of a slower decay of correlations in 
equilibrium. Their origin can be traced back to the subtle nonlinear effects described by mode-coupling 
theories. Measures of thermal conductivity with nonequilibrium molecular-dynamics method confirm the 
relevance of such effects for low-dimensional lattices even at very low temperatures. 



PACS. 63.10.+a Lattice dynamics: general theory 
duction 



1 Introduction 

The theoretical study of physical models in one spatial 
dimension is often justified by their mathematical sim- 
plicity jjj. Besides this, they display intriguing peculiar- 
ities that render them qualitatively different from their 
three-dimensional counterparts. Actually, a further and 
even more relevant motivation is the possibility of mod- 
ern experimental techniques to produce a variety of real 
systems that could , at least in principle, be effectively 
described by Id models. Furthermore, the latter are of 
great importance to approach the dynamical and statis- 
tical properties of biological molecules like DNA || or 
proteins ||. 

In the present paper we wish to study the problem of 
relaxation and heat transport in a classical chain of os- 
cillators with quadratic and cubic intersite potential. Al- 
though this is a truly textbook model of an insulating 
crystal, there is no, to the best of our knowledge, sys- 
tematic analysis of the effects of anharmonicity. The re- 
sults may be of interest to describe experimental systems 
like e.g. strongly anisotropic crystals @,||, polymers || 
or nanowires |j. Some theoretical investigations of ther- 
mal conductance for a quantum wire in ballistic [^) and 
anharmonic || regimes have been also recently presented. 

The first issue that we address is the following: which 
difference should one observe in this case with respect to 

a E-mail: lepri@avanzi.de.unifi.it 



05.60.-k Transport processes - 44.10.+i Heat con- 



the usual 3d systems? One may argue that the cubic non- 
linearity ( "three phonon" scattering) must be less efficient 
due to the stronger constraints to be fulfilled in Id, an 
argument that has been even invoked to explain several 
experimental results A signature of this fact on the 

dynamical correlations is thus to be expected. To sub- 
stantiate this argument, we studied the model within the 
framework of the memory- function formalism []To[ . A per- 
turbative calculation shows that correlation functions dis- 
play a power-law tail (Section II), that we study in some 
detail in the long-wavelength limit and compare, to some 
extent, with molecular dynamics simulations. The latter 
show that, even for relatively small anharmonicity, strong 
non-perturbative effects due to subtle nonlinear interac- 
tion arise for long enough chains. These can be quantita- 
tively accounted for by mode-coupling theories JlT| akin 
to the ones usually applied for dense fluids |l2] . 

The consequences of such a behaviour on nonequilib- 
rium energy transport is discussed in Section III. Even in 
absence of such strong memory effects, the thermal con- 
ductivity has been shown to diverge with the system size 
both for Id HOI and 2d lattices II- The character- 
istic divergence law, a signature of the mentioned mode- 
coupling effects Q , is recovered also in the present case. 
This signals again the overwhelming role of nonlinearity 
and fluctuations in low spatial dimensions. 
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2 Correlation functions for an anharmonic 
chain 

We consider a set of atoms of mass m arranged on a ring 
of N sites with spacing a. Let it; be the displacement of 
the Z-th particle from its equilibrium position la. The La- 
grangian reads as 
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(1) 

where g is the coupling constant. This is the lowest-order 
approximation of a generic anharmonic potential which 
couples nearest neighbours. For historical reasons, it is 
sometimes referred to as the Fermi-Pasta-Ulam a-model 
(FPU-a). For convenience, from now on we will always 
work in dimensionless units, where a, m and the angular 
frequency too are set to unity. This implies, for instance, 
that the sound speed atuo is also unity and that the energy 
and coupling constant (that we still denote with g) are 
measured respectively in units of mujoa 2 and mwofl • 

A difficulty of model (0) is the unboundedness of the 
potential. Therefore, in order to avoid runaway instabil- 
ity of trajectories (which could anyhow be overcome by 
adding even terms of higher order) we will deal with suf- 
ficiently small coupling constant and/or energies. 

Upon introducing the complex amplitudes U k through 
the discrete Fourier transform 
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where k is an integer ranging between — ^ + 1 and we 
can thus rewrite the equations of motion as 

U k + L0 2 k U k 
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(3) 

with Tk being the interaction force among modes. The 
condition on the indices of the sum is intended to be mod- 
ulo N while 

u> k = 2 | sin | . (4) 

is the usual harmonic (phononic) dispersion law. 



2.1 Memory effects 

Let us consider the normalized correlation function 

g k (t) = p<j*(u k (t)ui(o)) , (5) 

which is defined in such a way that G k {0) = 1 (/3 is the in- 
verse temperature). In the framework of the Mori-Zwanzig 



projection approach, it can be shown that it satisfies the 
equation of motion (Hi 



G k + r k (t - s)G k (s)ds + u>lG k 







(0) 



The memory function r k accounts for memory effects and 
can be connected to the nonlinear force by the fluctuation- 
dissipation relation 



r k (t) = fi(T k (t)T* k (0)) 



(J) 



Notice that here it is also implicitely assumed that slow 
terms possibly contained in T k are negligible in the ther- 
modynamic limit |p"i"| . 

In the following we will evaluate the memory function 
to the lowest order of perturbation theory. This amounts 
to evaluating the average in (|7|) on the measure of the un- 
perturbed system (g = 0). A straightforward calculation 
yields 



r k (t) 



Cg^ k 1 



N E 
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(8) 



where C is a suitable numerical constant. Rather remark- 
ably, the Laplace transform of r k can be evaluated exactly 
in the large N limit (see the Appendix for some details of 
the calculation). Indeed, if we let 2irk/N — > q and define 



r(q,z) 



we obtain the simple result 



r(q,t)e 



dt 



(9) 



r(q,z) = Kco 2 (q) 



(10) 

where we have introduced the new (small) coupling pa- 
rameter K = Cg 2 /2j3 and 
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Inversion of the Laplace transform yields |16| 
r(q,t) = Kto 2 (q) [J (f2 + (q)t) + J (f2-(q)t)) 



(11) 



(12) 



where Jo is the Bessel function. With this result at hand, 
we can give a closed expression for Q by solving Eq. (||) 
with the Laplace transform method that yields the formal 
solution (with Q(q, t = 0) = 0) 



iz + r(q, z) 



z 2 — uj 2 (q) — izT(q, z) 



(13) 



If, like in the present case, the dissipative effects are weak, 
the latter expression can be approximated as 



-i/2 



-i/2 



z - u(q) - ir{q, z)/2 z + uj{q) - ir(q, z)/2 

(14) 



Stefano Lepri: Memory Effects and Heat Transport in One-dimensional Insulators 



3 



and the correlation function are determined for every q by 
the standard methods of inverse Laplace transformation 
i.e. by determining the singularities in the complex plane 

of g. 



has been evaluated numerically finding that, as expected, 
the latter oscillates with a period close to 2tt while its 



envelope decays for large r as K r 



-3/2 



2.2 Long wavelengths 

Let us focus on the effective dynamics of the long-wavelength 
modes. The latter are the slowest ones and are directly 
responsible of energy transport and therefore of main in- 
terest here. For q — > one has from Eq. ([ll]) Q + w 4 
and u> w /2_ ~ M- This implies that the first term of the 
memory function decays on a much faster time scale and 
can be neglected so that 



r(q,z) 



Kq 2 



(15) 



Therefore, in this limit, the transform of the mode auto- 
correlation satisfies the scaling relation 



\q\ VM 



(16) 



where, from Eq. (|14| ) 
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(17) 



From Eqs.(|l6|) and (|T^) it follows immediately that G(q, t) 
depends only on the product \q\t. 

To estimate the behaviour of Q{q,t) in a more de- 
tailed way one has to study the properties of the func- 
tion G(w). Besides of the branch cut in z — ±1, the lat- 
ter has two simple poles in the upper part of the com- 
plex plane that, for small K, are approximatively given 
by ±1 + i^K 2 / 3 (we neglect the term 0(K 2 ^ 3 ) in the 
real part since this gives only a small correction to the 
unperturbed frequency) . Accordingly, the inversion of the 
Laplace transform ([ll]) yields two contributions to the de- 
cay, the exponential plus some power-law tails. More pre- 
cisely, the calculation yields 

g(q,t) cx exp(-7(flf)|t|) cosqt - Kf(K, \q\t) (18) 

where we have introduced the (short times) relaxation rate 

V3, 



7(9) 



i K 2/3 
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(19) 



Notice how the square-root in (y_7|) affects the functional 
dependence of the relaxation time on the coupling con- 
stant and therefore on temperature. The behaviour of the 
slowly decaying correction 



yl — w 2 
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(20) 



2.3 Comparison with the numerical simulations 

In the present section we compare the results of the per- 
turbative calculation with the outcomes of molecular dy- 
namics simulations. The latter were performed at equi- 
librium and in the microcanonical ensemble by integrat- 
ing the equations of motion with a third order symplec- 
tic algorithm flTf and starting from random initial condi- 
tions. We then let the system evolve for a certain transient 
time in order to start the measures from a generic phase- 
space point. The energy per particle has been fixed in ev- 
ery computation and the corresponding temperature has 
been measured as twice the average kinetic energy density. 
In computing spectra and correlations functions a Fast 
Fourier Transform routine has been used, and the data 
are usually averaged over an ensemble of several trajec- 
tories (typically between 20 and 200) to reduce statistical 
fluctuations. 

As a first check we compared the perturbative expres- 
sion of the memory function ( |l2| ) with the correlation 
of the nonlinear force appearing in Eq. (Q) (see Fig. 0). 
The two agree well on short times, but some deviation 
is observed after a few oscillation periods. Furthermore, 
the temperature and wavenumber dependence of the re- 
laxation rates has been checked by measuring the initial 
decay rate of the envelope of Q(q,t). As seen in Figs. || 
and |[ a reasonable agreement with the result of pertur- 
bation theory is obtained for q > 0.006 (corresponding 
to N < 10 3 ). However, for smaller wavenumbers (longer 
chains) the data crossover to a different behaviour that is 
consistent with the law q 5 / 3 . The latter is the result ex- 
pected from mode-coupling theory and has been numeri- 
cally observed for the model with quartic nonlinearity at 
much larger temperatures fli"|| . 

We can therefore conclude that, despite the smallness 
of the perturbative parameter (K <~ 0.01), the dynam- 
ical effects of nonlinear mode interaction take over on 
longer ( "hydrodynamic" ) time/length scales. Eventually, 
this originates substantial deviations from the results ob- 
tained in the previous section. 



3 Thermal conduction 

In the present Section we will study the consequences on 
the phenomenon of stationary heat transport along the 
chain. Let us discuss the case of linear response i.e. small 
thermal gradients. We first evaluate the asymptotic de- 
cay of the Green-Kubo integrand in the framework of the 
perturbative calculation described above and compare the 
results with numerical simulations. 
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Fig. 1. The normalized autocorrelation function of the nonlin- 
ear force Th for (3 = 10 4 , g = 1.0, N — 128, k = 1 (correspond- 
ing to q — 0.0490). The thick solid line is the (approximate) 
perturbative result Jo/2 (see eq.(|l2|)). 
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Fig. 3. Scaling of the relaxation rates with temperature for 
the first Fourier mode of a chain of length N = 256 (trian- 
gles) and N = 512 (squares). The straight line is the result of 
perturbation theory. 
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Fig. 2. Scaling of the relaxation rates with the wavenumber 
q for g = 0.25 and /3 = 10.0. The straight line is the result of 
perturbation theory. 

3.1 The Green-Kubo formula 

If we neglect the anharmonic term (which contributes to 
the conductivity only to order g 2 ), the heat current J is 
approximated by its harmonic part fl4|| 



e luJkt U~k- Taking into account the fact that ^ fc VkOJk\A^\ 2 = 
for symmetry reasons (vk — —V—k), one finds 

k 

The integrand appearing in the Green-Kubo formula for 
thermal conductivity 



k = ks/3 N / (J(t)J(0))dt 
Jo 



(23) 



contains correlation functions of quantities that are quadratic 
in Ak- As the latter quantities are slowly varying (Ak -C 
u>kAk) we get the approximate expression 

N{J(t)J(0)) « A-Y^vl{A k {t)Al{Q)) + c.c. (24) 

where we have also taken into account the equipartition 
of energy {uj 2 k \A k \ 2 ) = 1/(3. 

Let us assume that the large t behaviour of ( pi| ) is 
dominated by long- wavelength modes. On the basis of the 
above discussion one expects that in this limit the auto- 
correlation of A(q,t) will depend only on \q\t, i.e. 

(A(q,t)A*(q,0)) = ~A(\q\t) , (25) 



J 



Jh 



N/2 

— y 

k=-N/2+l 



where A is a suitable function. Replacing the sum in ( |24| ) 
VkUJk {UkUk ~ UkUkj (21) with an integral and taking into account the fact that 

v(q) —> 1 for q — > one gets 



where Vk = ui' k is the phase velocity. It is convenient to 
recast the above expression in the new variables Ak = 



N(J(t)J(0)) cc ± 



A{x)dx 



(26) 
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Since that the integral in this formula is convergent be- 
cause of the asymptotic behaviour of ( |20| ) , we can conclude 
that the Green-Kubo integrand decays as 1/t for large t 
and that the thermal conductivity is infinite. 

To estimate the divergence law for a finite chain, one 
can cut off the integral in the formula ( p3| ) up to some time 
proportional to the length N [Q. As a result, k should 
diverge logarithmically with the system size. In the next 
subsection we will compare this prediction with numerical 
results. 



3.2 Nonequilibrium simulations and finite-size 
conductivity 

A simple and efficient way to compute transport coeffi- 
cients is to use nonequilibrium molecular dynamics. First 
of all we consider the chain with fixed ends (uq = uyv+i = 
0) and let the first and the N-th oscillators interact with 
two reservoirs operating at different temperatures Tr — 
T+AT/2 and T L = T-AT/2, respectively. In such a way 
a net heat current flows through the lattice. The thermal 
conductivity for the finite chain can thus be measured as 
the ratio between such a flux and the applied temperature 
gradient. 

Among the several possible choices, we simulated the 
effect of the reservoirs by means of Nose-Hoover ther- 
mostatting method JlSf . The latter preserves the deter- 
ministic nature of the dynamics and is simply implemented 
by adding the force terms — (l U\ and —(run to the equa- 
tion of motion of the first and last oscillator respectively. 
The "thermal" variables (j,, Cr, evolve according to the 
dynamical equations Q 

c* = ^(f-i) , <m 

where r is the thermostat response time and controls the 
strength of the coupling between the reservoirs and the 
chain. The above prescriptions imply that the kinetic en- 
ergy of the boundary particles fluctuates around the im- 
posed average value, thus mimicking an interaction with 
a reservoir in canonical equilibrium. 

The simulations were performed integrating the equa- 
tion of motion for the bulk particles together with the ( p7| ) 
with a fourth-order Runge-Kutta algorithm. The bound- 
ary temperatures were kept fixed so that upon increasing 
the lattice length the applied temperature gradient AT/N 
is decreased and the Green-Kubo formula becomes more 
and more accurate. Both the time averaged kinetic tem- 



peratures Ti — iif and heat flux J where pj|,p0| 

J = + N+ 1 ~ Ul + -9(^+1 - u if] > 

(28) 

have been computed over a single trajectory (approxima- 
tively 10 6 time units) started from random initial condi- 
tions . In every run, a suitably long transient (about 10 4 
time units) has been discarded in order to let the system 
reach a statistically stationary state with each oscillator 
in local equilibrium. 

The thermal conductivity is thus computed as k — 
\J\N/AT. It has to be noted that the latter quantity rep- 
resents an effective transport coefficient including both 
boundary and bulk scattering mechanisms. Alternatively, 
one could as well consider, say, a subchain far enough from 
the ends and compute a bulk conductivity as the ratio be- 
tween | J\ and the actual temperature gradient measured 
there. Since the latter will also be inversely proportional 
to A, the resulting scaling with size (which is the main 
issue here) will be of course the same with both defini- 
tions. Besides this, the first choice avoids the difficulty of 
dealing with very small gradients (see below). 

The dependence of k on the chain length is reported in 
Fig. [| for two different series of simulations with T = 0.1 
and for relatively small (AT = 0.02) and large (AT = 0.1) 
applied temperature differences. In both cases n increases 
linearly with N for N < 10 3 . Moreover, no sizable tem- 
perature gradient forms along the chain. Both facts are a 
signature of the weakness of the anharmonic effects up to 
this time/length scales, as confirmed by comparing with 
data obtained (using the same setup and parameters) for 
the pure harmonic (g = 0) chain. The latter displays the 
expected linear increase of n with N |2lJ] and differ by less 
than a few percent from the anharmonic case in this range 
of system sizes (compare full dots and crosses in Fig. [|). 
The fact that k is smaller for larger AT can be thus be 
attributed to stronger boundary scattering that, in turn, 
tends to reduce the conductivity. 

For larger N the bulk anharmonic scattering takes over 
and, accordingly, the two curves for different AT approach 
one the other. The conductivity increases more slowly with 
the size, and we can tentatively measure the effective di- 
vergence law of the form N a . The data in the inset of 
Fig. H (full dots) are consistent with a convergence to the 
asympotic value a = 2/5 expected from the theory and 
observed in previous works on other Id models ]T^,p2[. 
Notice also that the crossover size is in reasonable agree- 
ment with the behavior reported in Fig. |[ 

In order to clarify the role of different anharmonic in- 
teractions, we also compared the above results with some 
measurements for a chain with quartic potential. More 
precisely, we considered the so-called Fermi-Pasta-Ulam 
(3— model (FPU-/?) where the cubic term in the lagrangian 
([l]) is replaced by (wz+i — it;) 4 /4. The corresponding cou- 
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pling constant has been set to the value 1.0 to have about 
the same ratio between the average anharmonic and har- 
monic potential energies as in the previous case (approx- 
imatively 0.04 for T = 0.1). With such a choice, the 
strength of the nonlinear terms (as measured by the suit- 
able perturbative parameters) are of the same order and 
a comparison between the two models makes sense. As 
shown again in Fig. |], the values of k are now definitely 
smaller (about one order of magnitude at N = 1024). 
Moreover, a linear temperature profile sets in along the 
chain and the data neatly approach a power-law behaviour 
with an exponent very close to 2/5 already for N > 64 (di- 
amonds in the inset of Fig. [I]). We therefore conclude that, 
consistently with the general argument given in the Intro- 
duction, the cubic nonlinearity is much less effective for 
what concerns the process of energy diffusion and trans- 
port in Id. 

In conclusion, we have shown that the dynamical cor- 
relations and transport in Id lattice are strongly affected 
by nonlinear effects which are not accounted for by the 
simple perturbative analysis reported above even in the 
weakly anharmonic regime. For the model with cubic in- 
tersite potential, they manifest themselves e.g. in a faster 
divergence (N a rather than logiV) of the thermal con- 
ductivity with the size N. Further consequences of those 
issues on other physical phenomena like energy diffusion or 
wavepacket propagation will be subject of future research. 

I acknowledge useful discussions with Roberto Livi, Anto- 
nio Politi, Stefano Ruffo and the research group Dynamics of 
Complex Systems in Florence. This work is partially supported 
by the INFM project Equilibrium and nonequilibrium dynamos 
in condensed matter. 
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Fig. 4. Thermal conductivity of the FPU-a model versus 
lattice length N for g = 0.25, T = 0.1, r = 1.0 and AT = 
0.1 (triangles), AT — 0.02 (full circles). The crosses (resp. 
diamonds) refer to the harmonic (resp. FPU-/3) models for T = 
0.1 and AT = 0.02. The inset shows the logarithmic derivative 
Alogn/AlogN versus logiV of the data for the FPU-a (full 
circles) and FPU-/3 (diamonds) respectively. 



It is convenient to perform the change of variable u = 
e »(<? -<?/2)/2 After some algebra, the integral reads as 



z 

2^ 



du 



Appendix: Evaluation of the memory function 



n 2 _(u + l/u) 2 /A + z 2 fl\ (u- l/u) 2 /A + z 2 

(31) 

where the frequencies fl± (q) are defined in Eq. ( pd| ) and the 
integration is along a path joining the points ± exp(— iq/2). 
By means of the further change of variable u 2 = Q, the in- 
tegration can be performed on the unit circle. The result 
( |l0[ ) is thus obtained by means of the theorem of residues. 



In this Appendix we sketch some details of the perturba- 
tive calculation of the memory function for the quadratic 
force. One wants to evaluate the Laplace transform of the 
sum appearing in Eq. (^). As a first step, we replace the 
sum for large N by the integral over the Brillouin zone 



1 



+7T 



— / dq' cosui(q')t cosu>(q — q')t 



(29) 



Its transform is thus given by 

C+7T 



1 



47T.L dq '\u{q')+u{q- q ')f-z 2 
1 



(w(g') - u{q - q')) 2 - z 2 



(30) 
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